In RodolfoPelinson/pelinson.et.al.2020: Data Analysis of Pelinson et al 2020
knitr::opts_chunk$set(echo = TRUE)
The goal of PredatorIsolationComm is to walk the user through the statistical analysis presented in:
"Pelinson et al 2020. Top predator introduction changes the effects of spatial isolation on freshwater community structure"
DOI: https://doi.org/10.1101/857318
You can install the last version of PredatorIsolationComm package from my GitHub with:
``` {r installing package not eval, eval = FALSE, echo = T}
install.packages("devtools")
devtools::install_github("RodolfoPelinson/PredatorIsolationComm")
library("PredatorIsolationComm")
This will give you access to all the data and functions used to produce the results shown in "Pelinson et al 2020. Top predator introduction changes the effects of spatial isolation on freshwater community structure".
``` {r installing package, eval = T, , echo = FALSE}
library("PredatorIsolationComm")
Other packages used here are:
vegan version 2.5-6
mvabund version 4.1.3
gllvm version 1.2.2
permute version 0.9-5
``` {r installing packages2, eval = T, results = "hide", warning = F, message = F}
library(vegan)
library(permute)
library(mvabund)
library(gllvm)
## Testing for differences in community structure
Loading necessary data.
```r
data(com,
com_SS1, fish_SS1,isolation_SS1, TRAITS_SS1,
com_SS2, fish_SS2,isolation_SS2, TRAITS_SS2,
com_SS3, fish_SS3,isolation_SS3, TRAITS_SS3)
When we considered all surveys together, we had to exclude ponds that were missing in the third survey from the second and first surveys to achieve a balanced design. The balanced design was necessary for our permutation procedure, which accounted by the non-independence of ponds sampled in different moments in time. The permutation procedure below shuffles rows freely ponds freely within the same pond IDs, and then pond IDs freely across treatments.
com_incomplete <- com[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4"),]
com_incomplete_oc <- decostand(com_incomplete, method = "pa")
com_incomplete <- com_incomplete[,which(colSums(com_incomplete_oc) > 3)]
isolation_incomplete <- isolation[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")]
fish_incomplete <- fish[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")]
survey_incomplete <- survey[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")]
ID_incomplete <- ID[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")]
ID_incomplete <- as.factor(as.character(ID_incomplete))
set.seed(3)
control <- permute::how(within = permute::Within(type = 'free'),
plots = Plots(strata = ID_incomplete, type = 'free'),
nperm = 10000)
permutations <- shuffleSet(nrow(com_incomplete), control = control)
Before running the models, we looked for the best probability distribution to model our community data (abundance). We considered Poisson and Negative Binomial distributions.
com_incomplete_mvabund <- mvabund(com_incomplete)
meanvar.plot(com_incomplete_mvabund, table =F, pch = 16)
abline(a = 0, b = 1, lwd = 2)
box(lwd = 2)
fit_incosistent_interaction_NB <- manyglm(com_incomplete_mvabund ~ survey_incomplete * (fish_incomplete * isolation_incomplete), family = "negative.binomial", cor.type = "I")
fit_incosistent_interaction_POIS <- manyglm(com_incomplete_mvabund ~ survey_incomplete * (fish_incomplete * isolation_incomplete), family = "poisson", cor.type = "I")
plot(fit_incosistent_interaction_POIS)
plot(fit_incosistent_interaction_NB)
We chose to use the Negative Binomial distribution to model data after inspecting plots of variances against means and the spread of Dunn-Smyth residuals against fitted values.
Now running the models and multivariate Likelihood Ratio Tests.
fit_no_effect <- manyglm(com_incomplete_mvabund ~ 1, family = "negative.binomial", cor.type = "I")
fit_Time <- manyglm(com_incomplete_mvabund ~ survey_incomplete, family = "negative.binomial", cor.type = "I")
fit_Time_fish <- manyglm(com_incomplete_mvabund ~ survey_incomplete + fish_incomplete, family = "negative.binomial", cor.type = "I")
fit_Time_fish_isolation <- manyglm(com_incomplete_mvabund ~ survey_incomplete + fish_incomplete + isolation_incomplete, family = "negative.binomial", cor.type = "I")
fit_Time_interaction <- manyglm(com_incomplete_mvabund ~ survey_incomplete + (fish_incomplete * isolation_incomplete), family = "negative.binomial", cor.type = "I")
fit_incosistent_interaction <- manyglm(com_incomplete_mvabund ~ survey_incomplete * (fish_incomplete * isolation_incomplete), family = "negative.binomial", cor.type = "I")
#Runing the Likelihood ratio tests
anova_incomplete <- anova(fit_no_effect,
fit_Time,
fit_Time_fish,
fit_Time_fish_isolation,
fit_Time_interaction,
fit_incosistent_interaction, bootID = permutations , test = "LR", resamp = "pit.trap")
anova_incomplete
Because we had significant three way interactions, we choose to run separete tests for each survey separetely.
First we have to change the abundance matrices into the right format fot mvabund and also prepare other matrices that will be used.
com_SS1_mvabund <- mvabund(com_SS1)
com_SS2_mvabund <- mvabund(com_SS2)
com_SS3_mvabund <- mvabund(com_SS3)
env_SS1 <- data.frame(fish_SS1, isolation_SS1)
env_SS2 <- data.frame(fish_SS2, isolation_SS2)
env_SS3 <- data.frame(fish_SS3, isolation_SS3)
com_oc <- decostand(com, method = "pa")
First Survey
fit_SS1_no_effect <- manyglm(com_SS1_mvabund ~ 1, family = "negative.binomial")
fit_SS1_fish <- manyglm(com_SS1_mvabund ~ fish_SS1, family = "negative.binomial")
fit_SS1_isolation <- manyglm(com_SS1_mvabund ~ isolation_SS1, family = "negative.binomial")
fit_SS1_interaction <- manyglm(com_SS1_mvabund ~ fish_SS1*isolation_SS1, family = "negative.binomial")
set.seed(3);anova_SS1 <- anova(fit_SS1_interaction, nBoot = 10000, p.uni = "none", test = "LR")
anova_SS1
There is a significant interaction between isolation and presence of fish.
Now, testing if the effect of treatments is mediated by trophic level:
env_SS1 <- data.frame(fish_SS1, isolation_SS1)
fit_SS1_no_trait_interaction <- traitglm(L = com_SS1, R = env_SS1, Q = TRAITS_SS1, formula = ~ isolation_SS1 * fish_SS1, method = "manyglm", col.intercepts = T)
fit_SS1_trait_pred_interaction <- traitglm(L = com_SS1, R = env_SS1, Q = TRAITS_SS1, formula = ~ (isolation_SS1:trophic) * (fish_SS1:trophic), method = "manyglm", col.intercepts = T)
set.seed(3);anova_SS1_trait_interaction <- anova.traitglm(fit_SS1_no_trait_interaction,
fit_SS1_trait_pred_interaction,
nBoot = 10000, test = "LR", show.time = "none")
anova_SS1_trait_interaction
No, it is not!
Plotting the ordination. We performed a model based unconstrained ordination using function gllvm from package gllvm and the plot_ordination function from this package.
fit_gllvm_no_effect_SS1<- gllvm(com_SS1,family = "negative.binomial", method = "VA", row.eff = F,n.init = 50, seed = 3, num.lv = 2)
colors_SS1 <- NULL
for (i in 1:dim(TRAITS_SS1)[1]){
if(TRAITS_SS1$trophic[i] == "Pr"){colors_SS1[i] <- "gold3"}
else {colors_SS1[i] <- "forestgreen"}
}
new_names_SS1 <- rownames(fit_gllvm_no_effect_SS1$params$theta)
for(i in 1:length(new_names_SS1)){
if(substring(new_names_SS1, 4,4)[i] == "_"){
new_names_SS1[i] <- substring(new_names_SS1[i], 5,7)
}
if(substring(new_names_SS1, 2,2)[i] == "_"){
new_names_SS1[i] <- substring(new_names_SS1[i], 3,5)
} else {new_names_SS1[i] <- substring(new_names_SS1[i], 1,3)}
}
plot_ordination(fit_gllvm_no_effect_SS1, x1 = fish_SS1, x2 = isolation_SS1,
species = T, elipse = T, sites = T,
site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"),
legend = T, legend_labels = c("Fishless - 30 m", "Fishless - 120 m","Fishless - 480 m","Fish - 30 m", "Fish - 120 m","Fish - 480 m"),
legend_ncol = 2, legend_horiz = F,species_col = colors_SS1, species_names = new_names_SS1, main = "First Survey", species_size = TRAITS_SS1$volume_log, name_species_size = TRAITS_SS1$volume_log*0.21)
Plotting pairwise comparisons coefficients and 95% confidence intervals for each species. There is a lot of code (I hope to improve it in the future!) to get the actual coefficients and their respective intervals. They are all in the file coefs_for_plot.R. We can call that file through the source function. You may have to download the file manually if you wish to run this code.
source("coefs_for_plots.R")
Now we can plot the coefficients using My_coefplot function from this package.
First, we compare the effect of fish in each isolation distance.
ab_SS1_pred <- order(colSums(com[which(TRAITS$trophic == "Pr")])[which(colSums(com_oc[which(survey == "1"),which(TRAITS$trophic == "Pr")]) > 2)], decreasing = T)
ab_SS1_cons <- order(colSums(com[which(TRAITS$trophic == "Non_Pred")])[which(colSums(com_oc[which(survey == "1"),which(TRAITS$trophic == "Non_Pred")]) > 2)], decreasing = T)
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(fish_effect_SS1_upper_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_upper_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(fish_effect_SS1_lower_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_lower_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "30 m - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
My_coefplot(c(as.matrix(c(fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(fish_effect_SS1_upper_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_upper_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(fish_effect_SS1_lower_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_lower_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "120 m - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
My_coefplot(c(as.matrix(c(fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(fish_effect_SS1_upper_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_upper_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(fish_effect_SS1_lower_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
fish_effect_SS1_lower_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "480 m - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
Now, the pairwise effect of isolation for ponds without fish.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "30 to 120 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
My_coefplot(c(as.matrix(c(isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "30 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
My_coefplot(c(as.matrix(c(isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "120 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
And, the pairwise effect of isolation for ponds with fish.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_fish_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_fish_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "30 to 120 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_fish_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_fish_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "30 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_fish_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
c(as.matrix(c(isolation_effect_fish_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred],
isolation_effect_fish_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))),
species_labels = c(names(isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]),
names(isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])),
xlab = "Effect on abundance",
main = "120 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
Second Survey
fit_SS2_no_effect <- manyglm(com_SS2_mvabund ~ 1, family = "negative.binomial")
fit_SS2_fish <- manyglm(com_SS2_mvabund ~ fish_SS2, family = "negative.binomial")
fit_SS2_fish_isolation <- manyglm(com_SS2_mvabund ~ fish_SS2+isolation_SS2, family = "negative.binomial")
fit_SS2_interaction <- manyglm(com_SS2_mvabund ~ fish_SS2*isolation_SS2, family = "negative.binomial")
set.seed(3);anova_SS2 <- anova(fit_SS2_interaction, nBoot = 10000, p.uni = "none", test = "LR")
anova_SS2
Again, we see a significant interaction between isolation and presence of fish.
Now, testing if the effect of treatments is mediated by trophic level:
fit_SS2_no_trait_interaction <- traitglm(L = com_SS2, R = env_SS2, Q = TRAITS_SS2, formula = ~ isolation_SS2 * fish_SS2, method = "manyglm", col.intercepts = T)
fit_SS2_trait_pred_interaction <- traitglm(L = com_SS2, R = env_SS2, Q = TRAITS_SS2, formula = ~ (isolation_SS2:trophic) * (fish_SS2:trophic), method = "manyglm", col.intercepts = T)
set.seed(3);anova_SS2_trait_interaction <- anova.traitglm(fit_SS2_no_trait_interaction,
fit_SS2_trait_pred_interaction,
nBoot = 10000, test = "LR", show.time = "none")
anova_SS2_trait_interaction
Yes, it is!
Plotting the ordination.
fit_gllvm_no_effect_SS2<- gllvm(com_SS2,family = "negative.binomial", method = "VA", row.eff = F,n.init = 50, seed = 3, num.lv = 2)
colors_SS2 <- NULL
for (i in 1:dim(TRAITS_SS2)[1]){
if(TRAITS_SS2$trophic[i] == "Pr"){colors_SS2[i] <- "gold3"}
else {colors_SS2[i] <- "forestgreen"}
}
new_names_SS2 <- rownames(fit_gllvm_no_effect_SS2$params$theta)
for(i in 1:length(new_names_SS2)){
if(substring(new_names_SS2, 4,4)[i] == "_"){
new_names_SS2[i] <- substring(new_names_SS2[i], 5,7)
}
if(substring(new_names_SS2, 2,2)[i] == "_"){
new_names_SS2[i] <- substring(new_names_SS2[i], 3,5)
} else {new_names_SS2[i] <- substring(new_names_SS2[i], 1,3)}
}
plot_ordination(fit_gllvm_no_effect_SS2, x1 = fish_SS2, x2 = isolation_SS2,
species = T, elipse = T, sites = T,
site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"),
legend = F, legend_labels = c("Fishless - 30 m", "Fishless - 120 m","Fishless - 480 m","Fish - 30 m", "Fish - 120 m","Fish - 480 m"),
legend_ncol = 2, legend_horiz = F,species_col = colors_SS2, species_names = new_names_SS2, main = "Second Survey", species_size = TRAITS_SS2$volume_log, name_species_size = TRAITS_SS2$volume_log*0.21)
Plotting coefficients and confidence intervals for each species.
First, we compare the effect of fish in each isolation distance.
ab_SS2_pred <- order(colSums(com[which(TRAITS$trophic == "Pr")])[which(colSums(com_oc[which(survey == "2"),which(TRAITS$trophic == "Pr")]) > 2)], decreasing = T)
ab_SS2_cons <- order(colSums(com[which(TRAITS$trophic == "Non_Pred")])[which(colSums(com_oc[which(survey == "2"),which(TRAITS$trophic == "Non_Pred")]) > 2)], decreasing = T)
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(fish_effect_SS2_upper_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_upper_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(fish_effect_SS2_lower_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_lower_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(fish_effect_SS2_upper_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_upper_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(fish_effect_SS2_lower_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_lower_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "120 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(fish_effect_SS2_upper_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_upper_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(fish_effect_SS2_lower_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
fish_effect_SS2_lower_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "480 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
Because the effect of treatments on aquatic insects was mediated by trophic level, we can assess the same parameters for predators and non-predators as a whole.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(rev(trait_fish_30_SS2))),
c(as.matrix(rev(trait_fish_30_SS2_upper))),
c(as.matrix(rev(trait_fish_30_SS2_lower))),
species_labels = c("Predators", "Non-Predators"),
xlab = NULL,
main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_fish_120_SS2))),
c(as.matrix(rev(trait_fish_120_SS2_upper))),
c(as.matrix(rev(trait_fish_120_SS2_lower))),
species_labels = c("Predators", "Non-Predators"),
xlab = NULL,
main = "120 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_fish_480_SS2))),
c(as.matrix(rev(trait_fish_480_SS2_upper))),
c(as.matrix(rev(trait_fish_480_SS2_lower))),
species_labels = c("Predators", "Non-Predators"),
xlab = NULL,
main = "480 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
Now, the pairwise effect of isolation for ponds without fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "30 to 120 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "30 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "120 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(rev(trait_30_120_SS2))),
c(as.matrix(rev(trait_30_120_SS2_upper))),
c(as.matrix(rev(trait_30_120_SS2_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 120 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_30_480_SS2))),
c(as.matrix(rev(trait_30_480_SS2_upper))),
c(as.matrix(rev(trait_30_480_SS2_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 480 m- SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_120_480_SS2))),
c(as.matrix(rev(trait_120_480_SS2_upper))),
c(as.matrix(rev(trait_120_480_SS2_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "120 m to 480 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
And, the pairwise effect of isolation for ponds with fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_fish_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_fish_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "30 to 120 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_fish_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_fish_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "30 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_fish_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
c(as.matrix(c(isolation_effect_fish_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred],
isolation_effect_fish_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))),
species_labels = c(names(isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]),
names(isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])),
xlab = "Effect on abundance",
main = "120 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(rev(trait_30_120_fish_SS2))),
c(as.matrix(rev(trait_30_120_fish_SS2_upper))),
c(as.matrix(rev(trait_30_120_fish_SS2_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 120 m - SS2 - FISH", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_30_480_fish_SS2))),
c(as.matrix(rev(trait_30_480_fish_SS2_upper))),
c(as.matrix(rev(trait_30_480_fish_SS2_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 480 m- SS2 - FISH", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_120_480_fish_SS2))),
c(as.matrix(rev(trait_120_480_fish_SS2_upper))),
c(as.matrix(rev(trait_120_480_fish_SS2_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "120 m to 480 m - SS2 - FISH", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
Third Survey
fit_SS3_no_effect <- manyglm(com_SS3_mvabund ~ 1, family = "negative.binomial")
fit_SS3_fish <- manyglm(com_SS3_mvabund ~ fish_SS3, family = "negative.binomial")
fit_SS3_fish_isolation <- manyglm(com_SS3_mvabund ~ fish_SS3+isolation_SS3, family = "negative.binomial")
fit_SS3_interaction <- manyglm(com_SS3_mvabund ~ fish_SS3*isolation_SS3, family = "negative.binomial")
set.seed(3);anova_SS3 <- anova(fit_SS3_interaction, nBoot = 10000, p.uni = "none", test = "LR")
anova_SS3
Again, we see a significant interaction between isolation and presence of fish.
Now, testing if the effect of treatments is mediated by trophic level:
fit_SS3_no_trait_interaction <- traitglm(L = com_SS3, R = env_SS3, formula = ~ isolation_SS3 * fish_SS3, method = "manyglm", col.intercepts = T)
fit_SS3_trait_pred_interaction <- traitglm(L = com_SS3, R = env_SS3, Q = TRAITS_SS3, formula = ~ (isolation_SS3:trophic) * (fish_SS3:trophic), method = "manyglm", col.intercepts = T)
set.seed(3);anova_SS3_trait_interaction <- anova(fit_SS3_no_trait_interaction,
fit_SS3_trait_pred_interaction,
nBoot = 10000, test = "LR", show.time = "none")
anova_SS3_trait_interaction
Yes! it is!
Plotting the ordination.
fit_gllvm_no_effect_SS3<- gllvm(com_SS3,family = "negative.binomial", method = "VA", row.eff = F,n.init = 50, seed = 3, num.lv = 2)
colors_SS3 <- NULL
for (i in 1:dim(TRAITS_SS3)[1]){
if(TRAITS_SS3$trophic[i] == "Pr"){colors_SS3[i] <- "gold3"}
else {colors_SS3[i] <- "forestgreen"}
}
new_names_SS3 <- rownames(fit_gllvm_no_effect_SS3$params$theta)
for(i in 1:length(new_names_SS3)){
if(substring(new_names_SS3, 4,4)[i] == "_"){
new_names_SS3[i] <- substring(new_names_SS3[i], 5,7)
}
if(substring(new_names_SS3, 2,2)[i] == "_"){
new_names_SS3[i] <- substring(new_names_SS3[i], 3,5)
} else {new_names_SS3[i] <- substring(new_names_SS3[i], 1,3)}
}
plot_ordination(fit_gllvm_no_effect_SS3, x1 = fish_SS3, x2 = isolation_SS3,
species = T, elipse = T, sites = T,
site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"),
legend = F, legend_labels = c("Fishless - 30 m", "Fishless - 120 m","Fishless - 480 m","Fish - 30 m", "Fish - 120 m","Fish - 480 m"),
legend_ncol = 2, legend_horiz = F,species_col = colors_SS3, species_names = new_names_SS3, main = "Third Survey", species_size = TRAITS_SS3$volume_log, name_species_size = TRAITS_SS3$volume_log*0.21, plot_centroid_dist = TRUE)
Plotting coefficients and confidence intervals for each species.
First, we compare the effect of fish in each isolation distance.
ab_SS3_pred <- order(colSums(com[which(TRAITS$trophic == "Pr")])[which(colSums(com_oc[which(survey == "3"),which(TRAITS$trophic == "Pr")]) > 2)], decreasing = T)
ab_SS3_cons <- order(colSums(com[which(TRAITS$trophic == "Non_Pred")])[which(colSums(com_oc[which(survey == "3"),which(TRAITS$trophic == "Non_Pred")]) > 2)], decreasing = T)
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(fish_effect_SS3_upper_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_upper_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(fish_effect_SS3_lower_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_lower_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(fish_effect_SS3_upper_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_upper_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(fish_effect_SS3_lower_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_lower_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "120 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(fish_effect_SS3_upper_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_upper_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(fish_effect_SS3_lower_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
fish_effect_SS3_lower_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "480 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
Because the effect of treatments on aquatic insects was mediated by trophic level, we can assess the same parameters for predators and non-predators as a whole.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(rev(trait_fish_30_SS3))),
c(as.matrix(rev(trait_fish_30_SS3_upper))),
c(as.matrix(rev(trait_fish_30_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_fish_120_SS3))),
c(as.matrix(rev(trait_fish_120_SS3_upper))),
c(as.matrix(rev(trait_fish_120_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_fish_480_SS3))),
c(as.matrix(rev(trait_fish_480_SS3_upper))),
c(as.matrix(rev(trait_fish_480_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
Now, the pairwise effect of isolation for ponds without fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "30 to 120 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "30 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "120 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(rev(trait_30_120_SS3))),
c(as.matrix(rev(trait_30_120_SS3_upper))),
c(as.matrix(rev(trait_30_120_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 120 m - SS3", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_30_480_SS3))),
c(as.matrix(rev(trait_30_480_SS3_upper))),
c(as.matrix(rev(trait_30_480_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 480 m- SS3", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_120_480_SS3))),
c(as.matrix(rev(trait_120_480_SS3_upper))),
c(as.matrix(rev(trait_120_480_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "120 m to 480 m - SS3", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
And, the pairwise effect of isolation for ponds with fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_fish_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_fish_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "30 to 120 - SS3", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_fish_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_fish_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "30 to 480 - SS3", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
My_coefplot(c(as.matrix(c(isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_fish_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
c(as.matrix(c(isolation_effect_fish_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred],
isolation_effect_fish_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))),
species_labels = c(names(isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]),
names(isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])),
xlab = "Effect on abundance",
main = "120 to 480 - SS3", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1)
My_coefplot(c(as.matrix(rev(trait_30_120_fish_SS3))),
c(as.matrix(rev(trait_30_120_fish_SS3_upper))),
c(as.matrix(rev(trait_30_120_fish_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 120 m - SS3 - FISH", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_30_480_fish_SS3))),
c(as.matrix(rev(trait_30_480_fish_SS3_upper))),
c(as.matrix(rev(trait_30_480_fish_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "30 m to 480 m- SS3 - FISH", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
My_coefplot(c(as.matrix(rev(trait_120_480_fish_SS3))),
c(as.matrix(rev(trait_120_480_fish_SS3_upper))),
c(as.matrix(rev(trait_120_480_fish_SS3_lower))),
species_labels = c("Predators", "Non Predators"),
xlab = NULL,
main = "120 m to 480 m - SS3 - FISH", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
Testing for differences in distances of centroids of fish and fishless treatments in different isolation treatments.
Here we will use a the function Dif_dist() to compute the distances between the centroids of fish and fishless treatments in each isolation treatment from the unconstrained ordination done by the gllvm function from package gllvm. Then, we compute the differences among the distances between the different distances and assess significance by permuting rows of the matrix always within fish and fishless treatments. By doing that we keep the mean differences among ponds with and without fish, but randomize any differences in the effect of fish across the isolation gradient.
First Survey
#Excluding samples for balanced design
set.seed(3)
Dif_dist_SS1<- Dif_dist(com = com_SS1, x1 = fish_SS1, x2 = isolation_SS1, nperm = 10000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 20, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F)
Dif_dist_SS1$distances
Dif_dist_SS1$diferences
Dif_dist_SS1$p.values
Second Survey
#Excluding samples for balanced design
set.seed(3)
Dif_dist_SS2<- Dif_dist(com = com_SS2, x1 = fish_SS2, x2 = isolation_SS1, nperm = 10000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 20, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F)
Dif_dist_SS2$distances
Dif_dist_SS2$diferences
Dif_dist_SS2$p.values
Third Survey
Here, again we had to exclude two ponds from the 30 and 120 isolation treatments to achieve a balanced design. We randomly selected one fishless pond from the 30m and 120m isolation treatments to be excluded.
set.seed(1)
first <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "30")], size = 1)
second <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "120")], size = 1)
com_SS3_incomplete <- com_SS3[-c(first,second),]
fish_SS3_incomplete <- fish_SS3[-c(first,second)]
isolation_SS3_incomplete <- isolation_SS3[-c(first,second)]
set.seed(1)
Dif_dist_SS3<- Dif_dist(com = com_SS3_incomplete, x1 = fish_SS3_incomplete,x2 = isolation_SS3_incomplete,nperm = 1000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 20, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F)
Dif_dist_SS3$distances
Dif_dist_SS3$diferences
Dif_dist_SS3$p.values
We can check if our permutation procedure is working by using the plot_null function to plot some of the null communities generated by Dif_dist(). This is an example of the first 15 permutations for the third survey:
set.seed(1)
par(mfrow = c(4,4))
plot_null(com_SS3_incomplete, nperm = 15 ,x1 = fish_SS3_incomplete, x2 = isolation_SS3_incomplete,
family = "negative.binomial", strata = T, orig_n.init = 10, xlim = c(-2.5,2.5), ylim = c(-2.5,2.5),
refit_perm = F, elipse = T, site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"))
It seems like everything is ok.
The p-values we got for the differences between centroids are valid for a specific combination of 3 out of 4 fishless ponds in 30 and 120 m distances. I had to do that to achieve a balanced design. We can check what is the proportion of times that the difference between fish and fishless ponds is significantly greater in more isolated ponds when compared to the lowest isolated (30m VS 480m) for all possible combinations. So we run the previous analysis again 1000 times (this is very time consuming) with different with different combinations of 3 out of 4 fishless ponds in 30 and 120 m distances.
#Excluding samples for balanced design
set.seed(1)
first <- rep(NA, 10)
second <- rep(NA, 10)
for(i in 1:1000){
first[i] <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "30")], size = 1)
second[i] <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "120")], size = 1)
}
dists_prop_SS3 <- array(NA, dim = c(6,3,1000))
dists_SS3 <- matrix(NA, nrow = 1000, ncol = 3)
for(i in 1:1000){
com_SS3_incomplete <- com_SS3[-c(first[i],second[i]),]
fish_SS3_incomplete <- fish_SS3[-c(first[i],second[i])]
isolation_SS3_incomplete <- isolation_SS3[-c(first[i],second[i])]
Dif_dist_SS3<- Dif_dist(com = com_SS3_incomplete, x1 = fish_SS3_incomplete,x2 = isolation_SS3_incomplete,nperm = 10000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 10, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F)
dists_prop_SS3[,,i] <- Dif_dist_SS3$p.values
dists_SS3[i,] <- c(as.matrix(Dif_dist_SS3$diferences))
message(paste("permutation ", i))
}
#Proportion of Greater
length(dists_prop_SS3[,3][dists_SS3[1,3,] > 0]) / length(dists_SS3[,3])
length(dists_prop_SS3[1,3,][dists_prop_SS3[1,3,] < 0.05]) / length(dists_prop_SS3[1,3,])
############################################
In 95% of cases the difference between fish and fishless ponds was greater in the 480m distance, if compared to the 30m distance. Also, this difference was significantly greater in isolated ponds for 75% of these different combinations.
RodolfoPelinson/pelinson.et.al.2020 documentation built on June 10, 2021, 5:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
The goal of PredatorIsolationComm is to walk the user through the statistical analysis presented in:
"Pelinson et al 2020. Top predator introduction changes the effects of spatial isolation on freshwater community structure"
DOI: https://doi.org/10.1101/857318
You can install the last version of PredatorIsolationComm package from my GitHub with:
``` {r installing package not eval, eval = FALSE, echo = T} install.packages("devtools") devtools::install_github("RodolfoPelinson/PredatorIsolationComm") library("PredatorIsolationComm")
This will give you access to all the data and functions used to produce the results shown in "Pelinson et al 2020. Top predator introduction changes the effects of spatial isolation on freshwater community structure".
``` {r installing package, eval = T, , echo = FALSE}
library("PredatorIsolationComm")
Other packages used here are:
vegan version 2.5-6
mvabund version 4.1.3
gllvm version 1.2.2
permute version 0.9-5
``` {r installing packages2, eval = T, results = "hide", warning = F, message = F} library(vegan) library(permute) library(mvabund) library(gllvm)
## Testing for differences in community structure
Loading necessary data.
```r
data(com,
com_SS1, fish_SS1,isolation_SS1, TRAITS_SS1,
com_SS2, fish_SS2,isolation_SS2, TRAITS_SS2,
com_SS3, fish_SS3,isolation_SS3, TRAITS_SS3)
When we considered all surveys together, we had to exclude ponds that were missing in the third survey from the second and first surveys to achieve a balanced design. The balanced design was necessary for our permutation procedure, which accounted by the non-independence of ponds sampled in different moments in time. The permutation procedure below shuffles rows freely ponds freely within the same pond IDs, and then pond IDs freely across treatments.
com_incomplete <- com[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4"),] com_incomplete_oc <- decostand(com_incomplete, method = "pa") com_incomplete <- com_incomplete[,which(colSums(com_incomplete_oc) > 3)] isolation_incomplete <- isolation[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")] fish_incomplete <- fish[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")] survey_incomplete <- survey[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")] ID_incomplete <- ID[which(ID != "A4" & ID != "B3" & ID != "C3" & ID != "C4")] ID_incomplete <- as.factor(as.character(ID_incomplete)) set.seed(3) control <- permute::how(within = permute::Within(type = 'free'), plots = Plots(strata = ID_incomplete, type = 'free'), nperm = 10000) permutations <- shuffleSet(nrow(com_incomplete), control = control)
Before running the models, we looked for the best probability distribution to model our community data (abundance). We considered Poisson and Negative Binomial distributions.
com_incomplete_mvabund <- mvabund(com_incomplete) meanvar.plot(com_incomplete_mvabund, table =F, pch = 16) abline(a = 0, b = 1, lwd = 2) box(lwd = 2) fit_incosistent_interaction_NB <- manyglm(com_incomplete_mvabund ~ survey_incomplete * (fish_incomplete * isolation_incomplete), family = "negative.binomial", cor.type = "I") fit_incosistent_interaction_POIS <- manyglm(com_incomplete_mvabund ~ survey_incomplete * (fish_incomplete * isolation_incomplete), family = "poisson", cor.type = "I") plot(fit_incosistent_interaction_POIS) plot(fit_incosistent_interaction_NB)
We chose to use the Negative Binomial distribution to model data after inspecting plots of variances against means and the spread of Dunn-Smyth residuals against fitted values.
Now running the models and multivariate Likelihood Ratio Tests.
fit_no_effect <- manyglm(com_incomplete_mvabund ~ 1, family = "negative.binomial", cor.type = "I") fit_Time <- manyglm(com_incomplete_mvabund ~ survey_incomplete, family = "negative.binomial", cor.type = "I") fit_Time_fish <- manyglm(com_incomplete_mvabund ~ survey_incomplete + fish_incomplete, family = "negative.binomial", cor.type = "I") fit_Time_fish_isolation <- manyglm(com_incomplete_mvabund ~ survey_incomplete + fish_incomplete + isolation_incomplete, family = "negative.binomial", cor.type = "I") fit_Time_interaction <- manyglm(com_incomplete_mvabund ~ survey_incomplete + (fish_incomplete * isolation_incomplete), family = "negative.binomial", cor.type = "I") fit_incosistent_interaction <- manyglm(com_incomplete_mvabund ~ survey_incomplete * (fish_incomplete * isolation_incomplete), family = "negative.binomial", cor.type = "I") #Runing the Likelihood ratio tests anova_incomplete <- anova(fit_no_effect, fit_Time, fit_Time_fish, fit_Time_fish_isolation, fit_Time_interaction, fit_incosistent_interaction, bootID = permutations , test = "LR", resamp = "pit.trap") anova_incomplete
Because we had significant three way interactions, we choose to run separete tests for each survey separetely.
First we have to change the abundance matrices into the right format fot mvabund and also prepare other matrices that will be used.
com_SS1_mvabund <- mvabund(com_SS1) com_SS2_mvabund <- mvabund(com_SS2) com_SS3_mvabund <- mvabund(com_SS3) env_SS1 <- data.frame(fish_SS1, isolation_SS1) env_SS2 <- data.frame(fish_SS2, isolation_SS2) env_SS3 <- data.frame(fish_SS3, isolation_SS3) com_oc <- decostand(com, method = "pa")
First Survey
fit_SS1_no_effect <- manyglm(com_SS1_mvabund ~ 1, family = "negative.binomial") fit_SS1_fish <- manyglm(com_SS1_mvabund ~ fish_SS1, family = "negative.binomial") fit_SS1_isolation <- manyglm(com_SS1_mvabund ~ isolation_SS1, family = "negative.binomial") fit_SS1_interaction <- manyglm(com_SS1_mvabund ~ fish_SS1*isolation_SS1, family = "negative.binomial") set.seed(3);anova_SS1 <- anova(fit_SS1_interaction, nBoot = 10000, p.uni = "none", test = "LR") anova_SS1
There is a significant interaction between isolation and presence of fish.
Now, testing if the effect of treatments is mediated by trophic level:
env_SS1 <- data.frame(fish_SS1, isolation_SS1) fit_SS1_no_trait_interaction <- traitglm(L = com_SS1, R = env_SS1, Q = TRAITS_SS1, formula = ~ isolation_SS1 * fish_SS1, method = "manyglm", col.intercepts = T) fit_SS1_trait_pred_interaction <- traitglm(L = com_SS1, R = env_SS1, Q = TRAITS_SS1, formula = ~ (isolation_SS1:trophic) * (fish_SS1:trophic), method = "manyglm", col.intercepts = T) set.seed(3);anova_SS1_trait_interaction <- anova.traitglm(fit_SS1_no_trait_interaction, fit_SS1_trait_pred_interaction, nBoot = 10000, test = "LR", show.time = "none") anova_SS1_trait_interaction
No, it is not!
Plotting the ordination. We performed a model based unconstrained ordination using function gllvm from package gllvm and the plot_ordination function from this package.
fit_gllvm_no_effect_SS1<- gllvm(com_SS1,family = "negative.binomial", method = "VA", row.eff = F,n.init = 50, seed = 3, num.lv = 2) colors_SS1 <- NULL for (i in 1:dim(TRAITS_SS1)[1]){ if(TRAITS_SS1$trophic[i] == "Pr"){colors_SS1[i] <- "gold3"} else {colors_SS1[i] <- "forestgreen"} } new_names_SS1 <- rownames(fit_gllvm_no_effect_SS1$params$theta) for(i in 1:length(new_names_SS1)){ if(substring(new_names_SS1, 4,4)[i] == "_"){ new_names_SS1[i] <- substring(new_names_SS1[i], 5,7) } if(substring(new_names_SS1, 2,2)[i] == "_"){ new_names_SS1[i] <- substring(new_names_SS1[i], 3,5) } else {new_names_SS1[i] <- substring(new_names_SS1[i], 1,3)} } plot_ordination(fit_gllvm_no_effect_SS1, x1 = fish_SS1, x2 = isolation_SS1, species = T, elipse = T, sites = T, site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"), legend = T, legend_labels = c("Fishless - 30 m", "Fishless - 120 m","Fishless - 480 m","Fish - 30 m", "Fish - 120 m","Fish - 480 m"), legend_ncol = 2, legend_horiz = F,species_col = colors_SS1, species_names = new_names_SS1, main = "First Survey", species_size = TRAITS_SS1$volume_log, name_species_size = TRAITS_SS1$volume_log*0.21)
Plotting pairwise comparisons coefficients and 95% confidence intervals for each species. There is a lot of code (I hope to improve it in the future!) to get the actual coefficients and their respective intervals. They are all in the file coefs_for_plot.R. We can call that file through the source function. You may have to download the file manually if you wish to run this code.
source("coefs_for_plots.R")
Now we can plot the coefficients using My_coefplot function from this package.
First, we compare the effect of fish in each isolation distance.
ab_SS1_pred <- order(colSums(com[which(TRAITS$trophic == "Pr")])[which(colSums(com_oc[which(survey == "1"),which(TRAITS$trophic == "Pr")]) > 2)], decreasing = T) ab_SS1_cons <- order(colSums(com[which(TRAITS$trophic == "Non_Pred")])[which(colSums(com_oc[which(survey == "1"),which(TRAITS$trophic == "Non_Pred")]) > 2)], decreasing = T) par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(fish_effect_SS1_upper_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_upper_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(fish_effect_SS1_lower_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_lower_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(fish_effect_SS1_30[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "30 m - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2) My_coefplot(c(as.matrix(c(fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(fish_effect_SS1_upper_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_upper_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(fish_effect_SS1_lower_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_lower_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(fish_effect_SS1_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "120 m - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2) My_coefplot(c(as.matrix(c(fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(fish_effect_SS1_upper_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_upper_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(fish_effect_SS1_lower_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], fish_effect_SS1_lower_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(fish_effect_SS1_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "480 m - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
Now, the pairwise effect of isolation for ponds without fish.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(isolation_effect_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "30 to 120 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2) My_coefplot(c(as.matrix(c(isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(isolation_effect_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "30 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2) My_coefplot(c(as.matrix(c(isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(isolation_effect_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "120 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
And, the pairwise effect of isolation for ponds with fish.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_fish_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_upper_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_fish_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_lower_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(isolation_effect_fish_SS1_30_120[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "30 to 120 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_fish_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_upper_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_fish_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_lower_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(isolation_effect_fish_SS1_30_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "30 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_fish_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_upper_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), c(as.matrix(c(isolation_effect_fish_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred], isolation_effect_fish_SS1_lower_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons]))), species_labels = c(names(isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Pr")][ab_SS1_pred]), names(isolation_effect_fish_SS1_120_480[which(TRAITS_SS1$trophic == "Non_Pred")][ab_SS1_cons])), xlab = "Effect on abundance", main = "120 to 480 - SS1", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 2)
Second Survey
fit_SS2_no_effect <- manyglm(com_SS2_mvabund ~ 1, family = "negative.binomial") fit_SS2_fish <- manyglm(com_SS2_mvabund ~ fish_SS2, family = "negative.binomial") fit_SS2_fish_isolation <- manyglm(com_SS2_mvabund ~ fish_SS2+isolation_SS2, family = "negative.binomial") fit_SS2_interaction <- manyglm(com_SS2_mvabund ~ fish_SS2*isolation_SS2, family = "negative.binomial") set.seed(3);anova_SS2 <- anova(fit_SS2_interaction, nBoot = 10000, p.uni = "none", test = "LR") anova_SS2
Again, we see a significant interaction between isolation and presence of fish.
Now, testing if the effect of treatments is mediated by trophic level:
fit_SS2_no_trait_interaction <- traitglm(L = com_SS2, R = env_SS2, Q = TRAITS_SS2, formula = ~ isolation_SS2 * fish_SS2, method = "manyglm", col.intercepts = T) fit_SS2_trait_pred_interaction <- traitglm(L = com_SS2, R = env_SS2, Q = TRAITS_SS2, formula = ~ (isolation_SS2:trophic) * (fish_SS2:trophic), method = "manyglm", col.intercepts = T) set.seed(3);anova_SS2_trait_interaction <- anova.traitglm(fit_SS2_no_trait_interaction, fit_SS2_trait_pred_interaction, nBoot = 10000, test = "LR", show.time = "none") anova_SS2_trait_interaction
Yes, it is!
Plotting the ordination.
fit_gllvm_no_effect_SS2<- gllvm(com_SS2,family = "negative.binomial", method = "VA", row.eff = F,n.init = 50, seed = 3, num.lv = 2) colors_SS2 <- NULL for (i in 1:dim(TRAITS_SS2)[1]){ if(TRAITS_SS2$trophic[i] == "Pr"){colors_SS2[i] <- "gold3"} else {colors_SS2[i] <- "forestgreen"} } new_names_SS2 <- rownames(fit_gllvm_no_effect_SS2$params$theta) for(i in 1:length(new_names_SS2)){ if(substring(new_names_SS2, 4,4)[i] == "_"){ new_names_SS2[i] <- substring(new_names_SS2[i], 5,7) } if(substring(new_names_SS2, 2,2)[i] == "_"){ new_names_SS2[i] <- substring(new_names_SS2[i], 3,5) } else {new_names_SS2[i] <- substring(new_names_SS2[i], 1,3)} } plot_ordination(fit_gllvm_no_effect_SS2, x1 = fish_SS2, x2 = isolation_SS2, species = T, elipse = T, sites = T, site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"), legend = F, legend_labels = c("Fishless - 30 m", "Fishless - 120 m","Fishless - 480 m","Fish - 30 m", "Fish - 120 m","Fish - 480 m"), legend_ncol = 2, legend_horiz = F,species_col = colors_SS2, species_names = new_names_SS2, main = "Second Survey", species_size = TRAITS_SS2$volume_log, name_species_size = TRAITS_SS2$volume_log*0.21)
Plotting coefficients and confidence intervals for each species.
First, we compare the effect of fish in each isolation distance.
ab_SS2_pred <- order(colSums(com[which(TRAITS$trophic == "Pr")])[which(colSums(com_oc[which(survey == "2"),which(TRAITS$trophic == "Pr")]) > 2)], decreasing = T) ab_SS2_cons <- order(colSums(com[which(TRAITS$trophic == "Non_Pred")])[which(colSums(com_oc[which(survey == "2"),which(TRAITS$trophic == "Non_Pred")]) > 2)], decreasing = T) par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(fish_effect_SS2_upper_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_upper_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(fish_effect_SS2_lower_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_lower_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(fish_effect_SS2_30[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(fish_effect_SS2_upper_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_upper_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(fish_effect_SS2_lower_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_lower_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(fish_effect_SS2_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "120 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(fish_effect_SS2_upper_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_upper_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(fish_effect_SS2_lower_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], fish_effect_SS2_lower_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(fish_effect_SS2_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "480 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
Because the effect of treatments on aquatic insects was mediated by trophic level, we can assess the same parameters for predators and non-predators as a whole.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(rev(trait_fish_30_SS2))), c(as.matrix(rev(trait_fish_30_SS2_upper))), c(as.matrix(rev(trait_fish_30_SS2_lower))), species_labels = c("Predators", "Non-Predators"), xlab = NULL, main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_fish_120_SS2))), c(as.matrix(rev(trait_fish_120_SS2_upper))), c(as.matrix(rev(trait_fish_120_SS2_lower))), species_labels = c("Predators", "Non-Predators"), xlab = NULL, main = "120 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_fish_480_SS2))), c(as.matrix(rev(trait_fish_480_SS2_upper))), c(as.matrix(rev(trait_fish_480_SS2_lower))), species_labels = c("Predators", "Non-Predators"), xlab = NULL, main = "480 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
Now, the pairwise effect of isolation for ponds without fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(isolation_effect_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "30 to 120 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(isolation_effect_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "30 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(isolation_effect_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "120 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(rev(trait_30_120_SS2))), c(as.matrix(rev(trait_30_120_SS2_upper))), c(as.matrix(rev(trait_30_120_SS2_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 120 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_30_480_SS2))), c(as.matrix(rev(trait_30_480_SS2_upper))), c(as.matrix(rev(trait_30_480_SS2_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 480 m- SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_120_480_SS2))), c(as.matrix(rev(trait_120_480_SS2_upper))), c(as.matrix(rev(trait_120_480_SS2_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "120 m to 480 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
And, the pairwise effect of isolation for ponds with fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_fish_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_upper_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_fish_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_lower_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(isolation_effect_fish_SS2_30_120[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "30 to 120 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_fish_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_upper_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_fish_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_lower_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(isolation_effect_fish_SS2_30_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "30 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_fish_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_upper_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), c(as.matrix(c(isolation_effect_fish_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred], isolation_effect_fish_SS2_lower_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons]))), species_labels = c(names(isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Pr")][ab_SS2_pred]), names(isolation_effect_fish_SS2_120_480[which(TRAITS_SS2$trophic == "Non_Pred")][ab_SS2_cons])), xlab = "Effect on abundance", main = "120 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(rev(trait_30_120_fish_SS2))), c(as.matrix(rev(trait_30_120_fish_SS2_upper))), c(as.matrix(rev(trait_30_120_fish_SS2_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 120 m - SS2 - FISH", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_30_480_fish_SS2))), c(as.matrix(rev(trait_30_480_fish_SS2_upper))), c(as.matrix(rev(trait_30_480_fish_SS2_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 480 m- SS2 - FISH", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_120_480_fish_SS2))), c(as.matrix(rev(trait_120_480_fish_SS2_upper))), c(as.matrix(rev(trait_120_480_fish_SS2_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "120 m to 480 m - SS2 - FISH", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
Third Survey
fit_SS3_no_effect <- manyglm(com_SS3_mvabund ~ 1, family = "negative.binomial") fit_SS3_fish <- manyglm(com_SS3_mvabund ~ fish_SS3, family = "negative.binomial") fit_SS3_fish_isolation <- manyglm(com_SS3_mvabund ~ fish_SS3+isolation_SS3, family = "negative.binomial") fit_SS3_interaction <- manyglm(com_SS3_mvabund ~ fish_SS3*isolation_SS3, family = "negative.binomial") set.seed(3);anova_SS3 <- anova(fit_SS3_interaction, nBoot = 10000, p.uni = "none", test = "LR") anova_SS3
Again, we see a significant interaction between isolation and presence of fish.
Now, testing if the effect of treatments is mediated by trophic level:
fit_SS3_no_trait_interaction <- traitglm(L = com_SS3, R = env_SS3, formula = ~ isolation_SS3 * fish_SS3, method = "manyglm", col.intercepts = T) fit_SS3_trait_pred_interaction <- traitglm(L = com_SS3, R = env_SS3, Q = TRAITS_SS3, formula = ~ (isolation_SS3:trophic) * (fish_SS3:trophic), method = "manyglm", col.intercepts = T) set.seed(3);anova_SS3_trait_interaction <- anova(fit_SS3_no_trait_interaction, fit_SS3_trait_pred_interaction, nBoot = 10000, test = "LR", show.time = "none") anova_SS3_trait_interaction
Yes! it is!
Plotting the ordination.
fit_gllvm_no_effect_SS3<- gllvm(com_SS3,family = "negative.binomial", method = "VA", row.eff = F,n.init = 50, seed = 3, num.lv = 2) colors_SS3 <- NULL for (i in 1:dim(TRAITS_SS3)[1]){ if(TRAITS_SS3$trophic[i] == "Pr"){colors_SS3[i] <- "gold3"} else {colors_SS3[i] <- "forestgreen"} } new_names_SS3 <- rownames(fit_gllvm_no_effect_SS3$params$theta) for(i in 1:length(new_names_SS3)){ if(substring(new_names_SS3, 4,4)[i] == "_"){ new_names_SS3[i] <- substring(new_names_SS3[i], 5,7) } if(substring(new_names_SS3, 2,2)[i] == "_"){ new_names_SS3[i] <- substring(new_names_SS3[i], 3,5) } else {new_names_SS3[i] <- substring(new_names_SS3[i], 1,3)} } plot_ordination(fit_gllvm_no_effect_SS3, x1 = fish_SS3, x2 = isolation_SS3, species = T, elipse = T, sites = T, site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"), legend = F, legend_labels = c("Fishless - 30 m", "Fishless - 120 m","Fishless - 480 m","Fish - 30 m", "Fish - 120 m","Fish - 480 m"), legend_ncol = 2, legend_horiz = F,species_col = colors_SS3, species_names = new_names_SS3, main = "Third Survey", species_size = TRAITS_SS3$volume_log, name_species_size = TRAITS_SS3$volume_log*0.21, plot_centroid_dist = TRUE)
Plotting coefficients and confidence intervals for each species.
First, we compare the effect of fish in each isolation distance.
ab_SS3_pred <- order(colSums(com[which(TRAITS$trophic == "Pr")])[which(colSums(com_oc[which(survey == "3"),which(TRAITS$trophic == "Pr")]) > 2)], decreasing = T) ab_SS3_cons <- order(colSums(com[which(TRAITS$trophic == "Non_Pred")])[which(colSums(com_oc[which(survey == "3"),which(TRAITS$trophic == "Non_Pred")]) > 2)], decreasing = T) par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(fish_effect_SS3_upper_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_upper_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(fish_effect_SS3_lower_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_lower_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(fish_effect_SS3_30[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(fish_effect_SS3_upper_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_upper_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(fish_effect_SS3_lower_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_lower_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(fish_effect_SS3_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "120 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(fish_effect_SS3_upper_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_upper_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(fish_effect_SS3_lower_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], fish_effect_SS3_lower_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(fish_effect_SS3_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "480 m - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
Because the effect of treatments on aquatic insects was mediated by trophic level, we can assess the same parameters for predators and non-predators as a whole.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(rev(trait_fish_30_SS3))), c(as.matrix(rev(trait_fish_30_SS3_upper))), c(as.matrix(rev(trait_fish_30_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_fish_120_SS3))), c(as.matrix(rev(trait_fish_120_SS3_upper))), c(as.matrix(rev(trait_fish_120_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_fish_480_SS3))), c(as.matrix(rev(trait_fish_480_SS3_upper))), c(as.matrix(rev(trait_fish_480_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m - SS2", cex.axis = 0.8, cex.main = 0.8,y_spa = 0.25, rect = T, rect_lim = 1)
Now, the pairwise effect of isolation for ponds without fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(isolation_effect_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "30 to 120 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(isolation_effect_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "30 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(isolation_effect_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "120 to 480 - SS2", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(rev(trait_30_120_SS3))), c(as.matrix(rev(trait_30_120_SS3_upper))), c(as.matrix(rev(trait_30_120_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 120 m - SS3", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_30_480_SS3))), c(as.matrix(rev(trait_30_480_SS3_upper))), c(as.matrix(rev(trait_30_480_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 480 m- SS3", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_120_480_SS3))), c(as.matrix(rev(trait_120_480_SS3_upper))), c(as.matrix(rev(trait_120_480_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "120 m to 480 m - SS3", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
And, the pairwise effect of isolation for ponds with fish for each species.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_fish_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_upper_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_fish_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_lower_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(isolation_effect_fish_SS3_30_120[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "30 to 120 - SS3", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_fish_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_upper_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_fish_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_lower_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(isolation_effect_fish_SS3_30_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "30 to 480 - SS3", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6) My_coefplot(c(as.matrix(c(isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_fish_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_upper_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), c(as.matrix(c(isolation_effect_fish_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred], isolation_effect_fish_SS3_lower_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons]))), species_labels = c(names(isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Pr")][ab_SS3_pred]), names(isolation_effect_fish_SS3_120_480[which(TRAITS_SS3$trophic == "Non_Pred")][ab_SS3_cons])), xlab = "Effect on abundance", main = "120 to 480 - SS3", cex.axis = 0.8, cex.main = 0.8, rect = T, rect_lim = 6)
For trophic level.
par(mfrow = c(1,3), mar = c(5, 5, 4, 2) + 0.1) My_coefplot(c(as.matrix(rev(trait_30_120_fish_SS3))), c(as.matrix(rev(trait_30_120_fish_SS3_upper))), c(as.matrix(rev(trait_30_120_fish_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 120 m - SS3 - FISH", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_30_480_fish_SS3))), c(as.matrix(rev(trait_30_480_fish_SS3_upper))), c(as.matrix(rev(trait_30_480_fish_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "30 m to 480 m- SS3 - FISH", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1) My_coefplot(c(as.matrix(rev(trait_120_480_fish_SS3))), c(as.matrix(rev(trait_120_480_fish_SS3_upper))), c(as.matrix(rev(trait_120_480_fish_SS3_lower))), species_labels = c("Predators", "Non Predators"), xlab = NULL, main = "120 m to 480 m - SS3 - FISH", cex.axis = 0.8, cex.main = 0.8, y_spa = 0.25, rect = T, rect_lim = 1)
Testing for differences in distances of centroids of fish and fishless treatments in different isolation treatments.
Here we will use a the function Dif_dist() to compute the distances between the centroids of fish and fishless treatments in each isolation treatment from the unconstrained ordination done by the gllvm function from package gllvm. Then, we compute the differences among the distances between the different distances and assess significance by permuting rows of the matrix always within fish and fishless treatments. By doing that we keep the mean differences among ponds with and without fish, but randomize any differences in the effect of fish across the isolation gradient.
First Survey
#Excluding samples for balanced design set.seed(3) Dif_dist_SS1<- Dif_dist(com = com_SS1, x1 = fish_SS1, x2 = isolation_SS1, nperm = 10000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 20, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F) Dif_dist_SS1$distances Dif_dist_SS1$diferences Dif_dist_SS1$p.values
Second Survey
#Excluding samples for balanced design set.seed(3) Dif_dist_SS2<- Dif_dist(com = com_SS2, x1 = fish_SS2, x2 = isolation_SS1, nperm = 10000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 20, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F) Dif_dist_SS2$distances Dif_dist_SS2$diferences Dif_dist_SS2$p.values
Third Survey
Here, again we had to exclude two ponds from the 30 and 120 isolation treatments to achieve a balanced design. We randomly selected one fishless pond from the 30m and 120m isolation treatments to be excluded.
set.seed(1) first <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "30")], size = 1) second <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "120")], size = 1) com_SS3_incomplete <- com_SS3[-c(first,second),] fish_SS3_incomplete <- fish_SS3[-c(first,second)] isolation_SS3_incomplete <- isolation_SS3[-c(first,second)] set.seed(1) Dif_dist_SS3<- Dif_dist(com = com_SS3_incomplete, x1 = fish_SS3_incomplete,x2 = isolation_SS3_incomplete,nperm = 1000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 20, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F) Dif_dist_SS3$distances Dif_dist_SS3$diferences Dif_dist_SS3$p.values
We can check if our permutation procedure is working by using the plot_null function to plot some of the null communities generated by Dif_dist(). This is an example of the first 15 permutations for the third survey:
set.seed(1) par(mfrow = c(4,4)) plot_null(com_SS3_incomplete, nperm = 15 ,x1 = fish_SS3_incomplete, x2 = isolation_SS3_incomplete, family = "negative.binomial", strata = T, orig_n.init = 10, xlim = c(-2.5,2.5), ylim = c(-2.5,2.5), refit_perm = F, elipse = T, site_colors = c("sienna1", "steelblue1","sienna3","steelblue3","sienna4", "steelblue4"))
It seems like everything is ok.
The p-values we got for the differences between centroids are valid for a specific combination of 3 out of 4 fishless ponds in 30 and 120 m distances. I had to do that to achieve a balanced design. We can check what is the proportion of times that the difference between fish and fishless ponds is significantly greater in more isolated ponds when compared to the lowest isolated (30m VS 480m) for all possible combinations. So we run the previous analysis again 1000 times (this is very time consuming) with different with different combinations of 3 out of 4 fishless ponds in 30 and 120 m distances.
#Excluding samples for balanced design set.seed(1) first <- rep(NA, 10) second <- rep(NA, 10) for(i in 1:1000){ first[i] <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "30")], size = 1) second[i] <- sample(row(com_SS3)[which(fish_SS3=="absent" & isolation_SS3 == "120")], size = 1) } dists_prop_SS3 <- array(NA, dim = c(6,3,1000)) dists_SS3 <- matrix(NA, nrow = 1000, ncol = 3) for(i in 1:1000){ com_SS3_incomplete <- com_SS3[-c(first[i],second[i]),] fish_SS3_incomplete <- fish_SS3[-c(first[i],second[i])] isolation_SS3_incomplete <- isolation_SS3[-c(first[i],second[i])] Dif_dist_SS3<- Dif_dist(com = com_SS3_incomplete, x1 = fish_SS3_incomplete,x2 = isolation_SS3_incomplete,nperm = 10000, family = "negative.binomial", num.lv = 2, strata = T, show.perm = F, orig_n.init = 10, perm_n.init = 1, type = "centroid", method = "VA", refit_perm = F) dists_prop_SS3[,,i] <- Dif_dist_SS3$p.values dists_SS3[i,] <- c(as.matrix(Dif_dist_SS3$diferences)) message(paste("permutation ", i)) }
#Proportion of Greater length(dists_prop_SS3[,3][dists_SS3[1,3,] > 0]) / length(dists_SS3[,3]) length(dists_prop_SS3[1,3,][dists_prop_SS3[1,3,] < 0.05]) / length(dists_prop_SS3[1,3,]) ############################################
In 95% of cases the difference between fish and fishless ponds was greater in the 480m distance, if compared to the 30m distance. Also, this difference was significantly greater in isolated ponds for 75% of these different combinations.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.